function r = hydro(Ti, TIME_STEPS)
  % First order hydro
  %  eliminating the u-dot-dot terms 
  %  fsolve search method
  close all;
  fig = figure; 
  load('cm', 'mycmap');
  set(fig, 'Colormap', mycmap);

  HYDROLIZE = false;

  GRID_SIZE = 16;
  DX = 0.01;
  Dt = 0.0001;
  VECTOR_SUBSAMPLING = 16;
  [gridx, gridy] = meshgrid([1:VECTOR_SUBSAMPLING:GRID_SIZE], [1:VECTOR_SUBSAMPLING:GRID_SIZE]);


  VIDEO = true;
  EVERY_FRAME=true;
  FPS = 30;
  SECONDS_PER_LIGHTYEAR = 15;
  if (VIDEO)
    vt = VideoWriter('v_t.avi', 'Motion JPEG AVI');
    vt.Quality = 100;
    open(vt);
  end

  ZERO_HYDRO = false;



  A = sqrt(3)*pi/18 - log(3)/2;

  T  = zeros(GRID_SIZE, GRID_SIZE);
  ux = zeros(GRID_SIZE, GRID_SIZE);
  uy = zeros(GRID_SIZE, GRID_SIZE);
  P1 = zeros(GRID_SIZE, GRID_SIZE);
  P2 = zeros(GRID_SIZE, GRID_SIZE);

  
  dtT  = zeros(GRID_SIZE, GRID_SIZE); 
  dtux = zeros(GRID_SIZE, GRID_SIZE); 
  dtuy = zeros(GRID_SIZE, GRID_SIZE); 
  dtP1 = zeros(GRID_SIZE, GRID_SIZE);
  dtP2 = zeros(GRID_SIZE, GRID_SIZE);

  Tdx = zeros(GRID_SIZE, GRID_SIZE); 
  Tdy = zeros(GRID_SIZE, GRID_SIZE); 
  uxdx = zeros(GRID_SIZE, GRID_SIZE); 
  uxdy = zeros(GRID_SIZE, GRID_SIZE); 
  uydx = zeros(GRID_SIZE, GRID_SIZE); 
  uydy = zeros(GRID_SIZE, GRID_SIZE); 
  P1dx = zeros(GRID_SIZE, GRID_SIZE); 
  P1dy = zeros(GRID_SIZE, GRID_SIZE); 
  P2dx = zeros(GRID_SIZE, GRID_SIZE); 
  P2dy = zeros(GRID_SIZE, GRID_SIZE); 


  % dxdxT = zeros(GRID_SIZE, GRID_SIZE); 
  % dxdyT = zeros(GRID_SIZE, GRID_SIZE); 
  % dydyT = zeros(GRID_SIZE, GRID_SIZE); 

  % dxdxux = zeros(GRID_SIZE, GRID_SIZE); 
  % dxdyux = zeros(GRID_SIZE, GRID_SIZE); 
  % dydyux = zeros(GRID_SIZE, GRID_SIZE); 

  % dxdxuy = zeros(GRID_SIZE, GRID_SIZE); 
  % dxdyuy = zeros(GRID_SIZE, GRID_SIZE); 
  % dydyuy = zeros(GRID_SIZE, GRID_SIZE); 



  sumT = 0;
  % initial conditions for vx vy come here
  % vx
  % vy
  % T
  TEMP_OFFSET =5;
  cx = 0.5 * GRID_SIZE + 0.5;
  cy = 0.5 * GRID_SIZE + 0.5;
  for x = 1:GRID_SIZE
    for y = 1:GRID_SIZE
       T(x, y) = exp(-((x-cx)^2 + (y-cy)^2)/4);
      % T(x,y) = 0.2*sin(2*pi*x/GRID_SIZE);
    end
  end
  % T(:,:) = Ti(:,:);
  T = T + TEMP_OFFSET;
  for x = 1:GRID_SIZE
     for y = 1:GRID_SIZE
       %T(x,y) = T(x,y) +  0.01 * rand();
       %vx(x,y) = 0.001 * rand();
       %vy(x,y) = 0.001 * rand();
     end
  end

  %vx = 0.000000 * ones(GRID_SIZE, GRID_SIZE);;
  %vy = 0.2 * ones(GRID_SIZE, GRID_SIZE);;
  
  % ****************************************************************
  %  hydrolize the initial conditions
  % ****************************************************************
  if (HYDROLIZE)
    level = 4;
    k = zeros(GRID_SIZE+1, GRID_SIZE+1);  
    [cx, cy] = meshgrid(0:GRID_SIZE-1,0:GRID_SIZE-1);
    [kx, ky] = meshgrid(-GRID_SIZE/2:GRID_SIZE/2);
    for kx = -GRID_SIZE/level:GRID_SIZE/level
      for ky = -GRID_SIZE/level:GRID_SIZE/level
        ft = T .* exp(-2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);
        k(kx+GRID_SIZE/2+1, ky+GRID_SIZE/2+1) = sum(sum(ft));
      end
    end
    T = zeros(GRID_SIZE:GRID_SIZE);
    for kx = -GRID_SIZE/level:GRID_SIZE/level
      for ky = -GRID_SIZE/level:GRID_SIZE/level
        T = T + k(kx+GRID_SIZE/2+1,ky+GRID_SIZE/2+1) * exp(2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);
        % vx = vx + 20*(rand()-0.5)*exp(2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);
        % vy = vy + 20*(rand()-0.5)*exp(2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);

      end
    end
    T = real(T / GRID_SIZE / GRID_SIZE);
    %ux = real(ux / GRID_SIZE / GRID_SIZE);
    %uy = real(uy / GRID_SIZE / GRID_SIZE);
  end



  % ****************************************************************
  %  compute derivative kernels
  % ****************************************************************
  % can upgrade on this -- trefethen pg 6

  d1 = zeros(GRID_SIZE, GRID_SIZE);
  d2 = zeros(GRID_SIZE, GRID_SIZE);
  for x = 1:GRID_SIZE
    for y=1:GRID_SIZE
      s1 = 0;
      s2 = 0;
      for k = -GRID_SIZE/2 : GRID_SIZE/2
        s1 = s1 + 2. * pi * i / GRID_SIZE / GRID_SIZE * k * exp(2 * pi * i * k * (x - y) / GRID_SIZE);
        s2 = s2 + (2. * pi * i / GRID_SIZE * k)^2 / GRID_SIZE * exp(2 * pi * i * k * (x - y) / GRID_SIZE);
      end
      d1(x,y) = real(s1)/DX;
      d2(x,y) = real(s2)/DX/DX;
    end
  end
  h = 2*pi/GRID_SIZE;
  col1 = [0 .5*(-1).^(1:GRID_SIZE-1).*cot((1:GRID_SIZE-1)*h/2)]  * h/DX;
  max(max(d1 - toeplitz(col1, col1([1 GRID_SIZE:-1:2]))))
  col2 = [(-pi^2/(3*h*h)-1/6) .5*(-1).^(0:GRID_SIZE-2).*(csc((1:GRID_SIZE-1)*h/2).^2)]  * h*h/(DX*DX);
  % d2 = toeplitz(col2);
  


  Tdx = dx(T);
  Tdy = dy(T);
  uxdx = dx(ux);
  uxdy = dy(ux);
  uydx = dx(uy);
  uydy = dy(uy);


  % **********************************************
  %  compute higher order initial conditions
  % **********************************************
%    Tdt = (-1).*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+uy.^2).^(-1).*(Tdx.*ux+Tdx.*ux.^3+T.*uxdx+Tdy.*uy+Tdy.*ux.^2.*uy+(-1).*T.*ux.*uxdy.*uy+Tdx.*ux.*uy.^2+T.*uxdx.*uy.^2+Tdy.*uy.^3+(-1).*T.*ux.*uy.*uydx+T.*uydy+T.*ux.^2.*uydy);
%    uxdt = T.^(-1).*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+uy.^2).^(-1).*((-2).*Tdx+(-2).*Tdx.*ux.^2+(-1).*T.*ux.*uxdx+(-1).*T.*ux.^3.*uxdx+(-1).*Tdy.*ux.*uy+(-2).*T.*uxdy.*uy+(-2).*T.*ux.^2.*uxdy.*uy+(-1).*Tdx.*uy.^2+(-1).*T.*uxdy.*uy.^3+(-1).*T.*ux.^2.*uy.*uydx+T.*ux.*uydy+T.*ux.^3.*uydy);
%    uydt = (-1).*T.^(-1).*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+uy.^2).^(-1).*(2.*Tdy+Tdy.*ux.^2+Tdx.*ux.*uy+(-1).*T.*uxdx.*uy+2.*Tdy.*uy.^2+T.*ux.*uxdy.*uy.^2+(-1).*T.*uxdx.*uy.^3+2.*T.*ux.*uydx+T.*ux.^3.*uydx+2.*T.*ux.*uy.^2.*uydx+T.*uy.*uydy+T.*uy.^3.*uydy);
%  
%  
%    P1 = (1/4).*T.^2.*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+uy.^2).^(-1).*((-2).*ux.*uxdt+(-1).*ux.^3.*uxdt+(-3).*ux.*uxdt.*uy.^2+(-2).*uxdx.*(1+ux.^2+uy.^2).^(1/2)+(-1).*ux.^2.*uxdx.*(1+ux.^2+uy.^2).^(1/2)+(-2).*ux.*uxdy.*uy.*(1+ux.^2+uy.^2).^(1/2)+(-1).*uxdx.*uy.^2.*(1+ux.^2+uy.^2).^(1/2)+2.*uy.*uydt+3.*ux.^2.*uy.*uydt+uy.^3.*uydt+2.*ux.*uy.*(1+ux.^2+uy.^2).^(1/2).*uydx+2.*(1+ux.^2+uy.^2).^(1/2).*uydy+ux.^2.*(1+ux.^2+uy.^2).^(1/2).*uydy+uy.^2.*(1+ux.^2+uy.^2).^(1/2).*uydy);
%  
%    P2 = (-1/2).*T.^2.*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+uy.^2).^(-1).*(uxdt.*uy+uxdt.*uy.^3+uxdy.*(1+ux.^2+uy.^2).^(1/2)+uxdy.*uy.^2.*(1+ux.^2+uy.^2).^(1/2)+ux.*uydt+ux.^3.*uydt+(1+ux.^2+uy.^2).^(1/2).*uydx+ux.^2.*(1+ux.^2+uy.^2).^(1/2).*uydx);




  P1dx = dx(P1);
  P1dy = dy(P1);
  P2dx = dx(P2);
  P2dy = dy(P2);

    t = 1;
    phydro();
    pause(0.003);

    %keyboard();
  % dxdxux = dxdx(ux);
  % dxdyux = dx(dy(ux));
  % dydyux = dydy(ux);
  % dxdxuy = dxdx(uy);
  % dxdyuy = dx(dy(uy));
  % dydyuy = dydy(uy);
  % dxdxT = dxdx(T);
  % dxdyT = dx(dy(T));
  % dydyT = dydy(T);



  dt = zeros(GRID_SIZE, GRID_SIZE*5);

  for t = 1:TIME_STEPS
    tic;

    % solving non-linear equation of type:
    % L(u, T, u-dot, T-dot) + R(u, T) == 0  -- 3 components
    % each of L, R has 3 components 

    % precompute the R functions
    % R1
     
    R1 = (1/2).*(1+ux.^2+uy.^2).^(-3/2).*(4.*P1dx.*ux+4.*P2dy.*ux+3.*T.^2.*Tdx.*ux+4.*P1dx.*ux.^3+6.*P2dy.*ux.^3+3.*T.^2.*Tdx.*ux.^3+2.*P2dy.*ux.^5+4.*P1.*uxdx+3.*T.^3.*uxdx+4.*P1.*ux.^2.*uxdx+3.*T.^3.*ux.^2.*uxdx+4.*P2.*uxdy+6.*P2.*ux.^2.*uxdy+4.*P2.*ux.^4.*uxdy+(-4).*P1dy.*uy+4.*P2dx.*uy+3.*T.^2.*Tdy.*uy+(-8).*P1dy.*ux.^2.*uy+2.*P2dx.*ux.^2.*uy+3.*T.^2.*Tdy.*ux.^2.*uy+(-4).*P1dy.*ux.^4.*uy+(-2).*P2dx.*ux.^4.*uy+(-2).*P2.*ux.^3.*uxdx.*uy+(-4).*P1.*ux.*uxdy.*uy+(-3).*T.^3.*ux.*uxdy.*uy+(-8).*P1.*ux.^3.*uxdy.*uy+(-3).*T.^3.*ux.^3.*uxdy.*uy+8.*P1dx.*ux.*uy.^2+2.*P2dy.*ux.*uy.^2+3.*T.^2.*Tdx.*ux.*uy.^2+4.*P1dx.*ux.^3.*uy.^2+4.*P1.*uxdx.*uy.^2+6.*T.^3.*uxdx.*uy.^2+4.*P1.*ux.^2.*uxdx.*uy.^2+3.*T.^3.*ux.^2.*uxdx.*uy.^2+2.*P2.*uxdy.*uy.^2+(-2).*P2.*ux.^2.*uxdy.*uy.^2+(-4).*P1dy.*uy.^3+6.*P2dx.*uy.^3+3.*T.^2.*Tdy.*uy.^3+(-4).*P1dy.*ux.^2.*uy.^3+2.*P2.*ux.*uxdx.*uy.^3+(-4).*P1.*ux.*uxdy.*uy.^3+(-3).*T.^3.*ux.*uxdy.*uy.^3+4.*P1dx.*ux.*uy.^4+(-2).*P2dy.*ux.*uy.^4+3.*T.^3.*uxdx.*uy.^4+(-2).*P2.*uxdy.*uy.^4+2.*P2dx.*uy.^5+4.*P2.*uydx+2.*P2.*ux.^2.*uydx+(-2).*P2.*ux.^4.*uydx+4.*P1.*ux.*uy.*uydx+(-3).*T.^3.*ux.*uy.*uydx+4.*P1.*ux.^3.*uy.*uydx+(-3).*T.^3.*ux.^3.*uy.*uydx+6.*P2.*uy.^2.*uydx+(-2).*P2.*ux.^2.*uy.^2.*uydx+8.*P1.*ux.*uy.^3.*uydx+(-3).*T.^3.*ux.*uy.^3.*uydx+4.*P2.*uy.^4.*uydx+(-4).*P1.*uydy+3.*T.^3.*uydy+(-4).*P1.*ux.^2.*uydy+6.*T.^3.*ux.^2.*uydy+3.*T.^3.*ux.^4.*uydy+2.*P2.*ux.^3.*uy.*uydy+(-4).*P1.*uy.^2.*uydy+3.*T.^3.*uy.^2.*uydy+(-4).*P1.*ux.^2.*uy.^2.*uydy+3.*T.^3.*ux.^2.*uy.^2.*uydy+(-2).*P2.*ux.*uy.^3.*uydy);

    R2 = (1/2).*(1+ux.^2+uy.^2).^(-2).*(4.*P1dx+4.*P2dy+3.*T.^2.*Tdx+8.*P1dx.*ux.^2+10.*P2dy.*ux.^2+6.*T.^2.*Tdx.*ux.^2+4.*P1dx.*ux.^4+8.*P2dy.*ux.^4+3.*T.^2.*Tdx.*ux.^4+2.*P2dy.*ux.^6+4.*P1.*ux.*uxdx+3.*T.^3.*ux.*uxdx+4.*P1.*ux.^3.*uxdx+3.*T.^3.*ux.^3.*uxdx+4.*P2.*ux.*uxdy+8.*P2.*ux.^3.*uxdy+4.*P2.*ux.^5.*uxdy+(-4).*P1dy.*ux.*uy+(-8).*P1dy.*ux.^3.*uy+(-2).*P2dx.*ux.^3.*uy+(-4).*P1dy.*ux.^5.*uy+(-2).*P2dx.*ux.^5.*uy+4.*P2.*uxdx.*uy+(-2).*P2.*ux.^4.*uxdx.*uy+(-4).*P1.*uxdy.*uy+(-12).*P1.*ux.^2.*uxdy.*uy+(-3).*T.^3.*ux.^2.*uxdy.*uy+(-8).*P1.*ux.^4.*uxdy.*uy+(-3).*T.^3.*ux.^4.*uxdy.*uy+8.*P1dx.*uy.^2+6.*P2dy.*uy.^2+6.*T.^2.*Tdx.*uy.^2+12.*P1dx.*ux.^2.*uy.^2+6.*P2dy.*ux.^2.*uy.^2+6.*T.^2.*Tdx.*ux.^2.*uy.^2+4.*P1dx.*ux.^4.*uy.^2+12.*P1.*ux.*uxdx.*uy.^2+6.*T.^3.*ux.*uxdx.*uy.^2+4.*P1.*ux.^3.*uxdx.*uy.^2+3.*T.^3.*ux.^3.*uxdx.*uy.^2+2.*P2.*ux.^3.*uxdy.*uy.^2+(-4).*P1dy.*ux.*uy.^3+2.*P2dx.*ux.*uy.^3+(-4).*P1dy.*ux.^3.*uy.^3+8.*P2.*uxdx.*uy.^3+(-2).*P2.*ux.^2.*uxdx.*uy.^3+(-4).*P1.*uxdy.*uy.^3+(-12).*P1.*ux.^2.*uxdy.*uy.^3+(-3).*T.^3.*ux.^2.*uxdy.*uy.^3+4.*P1dx.*uy.^4+2.*P2dy.*uy.^4+3.*T.^2.*Tdx.*uy.^4+4.*P1dx.*ux.^2.*uy.^4+(-2).*P2dy.*ux.^2.*uy.^4+8.*P1.*ux.*uxdx.*uy.^4+3.*T.^3.*ux.*uxdx.*uy.^4+(-6).*P2.*ux.*uxdy.*uy.^4+2.*P2dx.*ux.*uy.^5+4.*P2.*uxdx.*uy.^5+(-2).*P2.*ux.^3.*uydx+(-2).*P2.*ux.^5.*uydx+4.*P1.*ux.^2.*uy.*uydx+(-3).*T.^3.*ux.^2.*uy.*uydx+4.*P1.*ux.^4.*uy.*uydx+(-3).*T.^3.*ux.^4.*uy.*uydx+2.*P2.*ux.*uy.^2.*uydx+2.*P2.*ux.^3.*uy.^2.*uydx+(-3).*T.^3.*ux.^2.*uy.^3.*uydx+3.*T.^3.*ux.*uydy+6.*T.^3.*ux.^3.*uydy+3.*T.^3.*ux.^5.*uydy+(-2).*P2.*ux.^2.*uy.*uydy+(-2).*P2.*ux.^4.*uy.*uydy+4.*P1.*ux.*uy.^2.*uydy+3.*T.^3.*ux.*uy.^2.*uydy+4.*P1.*ux.^3.*uy.^2.*uydy+3.*T.^3.*ux.^3.*uy.^2.*uydy+2.*P2.*uy.^3.*uydy+2.*P2.*ux.^2.*uy.^3.*uydy);

    R3 = (1/2).*(1+ux.^2+uy.^2).^(-2).*((-4).*P1dy+4.*P2dx+3.*T.^2.*Tdy+(-8).*P1dy.*ux.^2+6.*P2dx.*ux.^2+6.*T.^2.*Tdy.*ux.^2+(-4).*P1dy.*ux.^4+2.*P2dx.*ux.^4+3.*T.^2.*Tdy.*ux.^4+2.*P2.*ux.^3.*uxdx+4.*P1dx.*ux.*uy+4.*P1dx.*ux.^3.*uy+2.*P2dy.*ux.^3.*uy+2.*P2dy.*ux.^5.*uy+3.*T.^3.*uxdx.*uy+(-4).*P1.*ux.^2.*uxdx.*uy+3.*T.^3.*ux.^2.*uxdx.*uy+2.*P2.*ux.^2.*uxdy.*uy+(-8).*P1dy.*uy.^2+10.*P2dx.*uy.^2+6.*T.^2.*Tdy.*uy.^2+(-12).*P1dy.*ux.^2.*uy.^2+6.*P2dx.*ux.^2.*uy.^2+6.*T.^2.*Tdy.*ux.^2.*uy.^2+(-4).*P1dy.*ux.^4.*uy.^2+(-2).*P2dx.*ux.^4.*uy.^2+(-2).*P2.*ux.*uxdx.*uy.^2+2.*P2.*ux.^3.*uxdx.*uy.^2+(-4).*P1.*ux.*uxdy.*uy.^2+(-3).*T.^3.*ux.*uxdy.*uy.^2+(-3).*T.^3.*ux.^3.*uxdy.*uy.^2+8.*P1dx.*ux.*uy.^3+(-2).*P2dy.*ux.*uy.^3+4.*P1dx.*ux.^3.*uy.^3+6.*T.^3.*uxdx.*uy.^3+(-4).*P1.*ux.^2.*uxdx.*uy.^3+3.*T.^3.*ux.^2.*uxdx.*uy.^3+(-2).*P2.*uxdy.*uy.^3+2.*P2.*ux.^2.*uxdy.*uy.^3+(-4).*P1dy.*uy.^4+8.*P2dx.*uy.^4+3.*T.^2.*Tdy.*uy.^4+(-4).*P1dy.*ux.^2.*uy.^4+(-2).*P2.*ux.*uxdx.*uy.^4+(-4).*P1.*ux.*uxdy.*uy.^4+(-3).*T.^3.*ux.*uxdy.*uy.^4+4.*P1dx.*ux.*uy.^5+(-2).*P2dy.*ux.*uy.^5+3.*T.^3.*uxdx.*uy.^5+(-2).*P2.*uxdy.*uy.^5+2.*P2dx.*uy.^6+4.*P1.*ux.*uydx+4.*P1.*ux.^3.*uydx+4.*P2.*uy.*uydx+(-6).*P2.*ux.^4.*uy.*uydx+12.*P1.*ux.*uy.^2.*uydx+(-3).*T.^3.*ux.*uy.^2.*uydx+12.*P1.*ux.^3.*uy.^2.*uydx+(-3).*T.^3.*ux.^3.*uy.^2.*uydx+8.*P2.*uy.^3.*uydx+2.*P2.*ux.^2.*uy.^3.*uydx+8.*P1.*ux.*uy.^4.*uydx+(-3).*T.^3.*ux.*uy.^4.*uydx+4.*P2.*uy.^5.*uydx+4.*P2.*ux.*uydy+8.*P2.*ux.^3.*uydy+4.*P2.*ux.^5.*uydy+(-4).*P1.*uy.*uydy+3.*T.^3.*uy.*uydy+(-12).*P1.*ux.^2.*uy.*uydy+6.*T.^3.*ux.^2.*uy.*uydy+(-8).*P1.*ux.^4.*uy.*uydy+3.*T.^3.*ux.^4.*uy.*uydy+(-2).*P2.*ux.^3.*uy.^2.*uydy+(-4).*P1.*uy.^3.*uydy+3.*T.^3.*uy.^3.*uydy+(-4).*P1.*ux.^2.*uy.^3.*uydy+3.*T.^3.*ux.^2.*uy.^3.*uydy+(-2).*P2.*ux.*uy.^4.*uydy);

    R4 = (-1/2).*T.^(-1).*(1+ux.^2+uy.^2).^(-1).*((-4).*P1.*T+(-8).*P1.*T.*ux.^2+(-4).*P1.*T.*ux.^4+(-2).*P1.*uxdx+(-2).*A.*P1.*uxdx+(-1).*T.^3.*uxdx+(-2).*P1.*ux.^2.*uxdx+(-2).*A.*P1.*ux.^2.*uxdx+(-1).*T.^3.*ux.^2.*uxdx+4.*A.*P2.*uxdy+6.*A.*P2.*ux.^2.*uxdy+2.*A.*P2.*ux.^4.*uxdy+(-4).*P2.*T.*ux.*uy+(-4).*P2.*T.*ux.^3.*uy+(-2).*P2.*ux.*uxdx.*uy+(-2).*A.*P2.*ux.*uxdx.*uy+2.*P1.*ux.*uxdy.*uy+(-2).*A.*P1.*ux.*uxdy.*uy+(-1).*T.^3.*ux.*uxdy.*uy+2.*P1.*ux.^3.*uxdy.*uy+(-2).*A.*P1.*ux.^3.*uxdy.*uy+(-1).*T.^3.*ux.^3.*uxdy.*uy+(-4).*P1.*T.*uy.^2+(-4).*P1.*T.*ux.^2.*uy.^2+(-2).*P1.*uxdx.*uy.^2+(-2).*A.*P1.*uxdx.*uy.^2+(-1).*T.^3.*uxdx.*uy.^2+(-2).*P1.*ux.^2.*uxdx.*uy.^2+(-2).*A.*P1.*ux.^2.*uxdx.*uy.^2+T.^3.*ux.^2.*uxdx.*uy.^2+2.*A.*P2.*uxdy.*uy.^2+2.*P2.*ux.^2.*uxdy.*uy.^2+(-4).*P2.*T.*ux.*uy.^3+(-2).*P2.*ux.*uxdx.*uy.^3+(-2).*A.*P2.*ux.*uxdx.*uy.^3+(-4).*A.*P2.*uydx+(-6).*A.*P2.*ux.^2.*uydx+(-2).*A.*P2.*ux.^4.*uydx+2.*P1.*ux.*uy.*uydx+6.*A.*P1.*ux.*uy.*uydx+(-1).*T.^3.*ux.*uy.*uydx+2.*P1.*ux.^3.*uy.*uydx+6.*A.*P1.*ux.^3.*uy.*uydx+(-1).*T.^3.*ux.^3.*uy.*uydx+(-2).*A.*P2.*uy.^2.*uydx+2.*P2.*ux.^2.*uy.^2.*uydx+4.*A.*P2.*ux.^2.*uy.^2.*uydx+(-2).*P1.*uydy+(-2).*A.*P1.*uydy+T.^3.*uydy+(-4).*P1.*ux.^2.*uydy+(-4).*A.*P1.*ux.^2.*uydy+2.*T.^3.*ux.^2.*uydy+(-2).*P1.*ux.^4.*uydy+(-2).*A.*P1.*ux.^4.*uydy+T.^3.*ux.^4.*uydy+(-2).*P2.*ux.*uy.*uydy+(-2).*A.*P2.*ux.*uy.*uydy+(-2).*P2.*ux.^3.*uy.*uydy+(-2).*A.*P2.*ux.^3.*uy.*uydy);

    R5 = (-1/2).*T.^(-1).*(1+ux.^2+uy.^2).^(-1).*(4.*P1.*T+4.*P1.*T.*ux.^2+2.*P1.*uxdx+2.*A.*P1.*uxdx+T.^3.*uxdx+(-4).*A.*P2.*uxdy+(-2).*A.*P2.*ux.^2.*uxdy+(-4).*P2.*T.*ux.*uy+(-4).*P2.*T.*ux.^3.*uy+(-2).*P2.*ux.*uxdx.*uy+(-2).*A.*P2.*ux.*uxdx.*uy+(-2).*P1.*ux.*uxdy.*uy+(-6).*A.*P1.*ux.*uxdy.*uy+(-1).*T.^3.*ux.*uxdy.*uy+8.*P1.*T.*uy.^2+4.*P1.*T.*ux.^2.*uy.^2+4.*P1.*uxdx.*uy.^2+4.*A.*P1.*uxdx.*uy.^2+2.*T.^3.*uxdx.*uy.^2+(-6).*A.*P2.*uxdy.*uy.^2+2.*P2.*ux.^2.*uxdy.*uy.^2+4.*A.*P2.*ux.^2.*uxdy.*uy.^2+(-4).*P2.*T.*ux.*uy.^3+(-2).*P2.*ux.*uxdx.*uy.^3+(-2).*A.*P2.*ux.*uxdx.*uy.^3+(-2).*P1.*ux.*uxdy.*uy.^3+(-6).*A.*P1.*ux.*uxdy.*uy.^3+(-1).*T.^3.*ux.*uxdy.*uy.^3+4.*P1.*T.*uy.^4+2.*P1.*uxdx.*uy.^4+2.*A.*P1.*uxdx.*uy.^4+T.^3.*uxdx.*uy.^4+(-2).*A.*P2.*uxdy.*uy.^4+4.*A.*P2.*uydx+2.*A.*P2.*ux.^2.*uydx+(-2).*P1.*ux.*uy.*uydx+2.*A.*P1.*ux.*uy.*uydx+(-1).*T.^3.*ux.*uy.*uydx+6.*A.*P2.*uy.^2.*uydx+2.*P2.*ux.^2.*uy.^2.*uydx+(-2).*P1.*ux.*uy.^3.*uydx+2.*A.*P1.*ux.*uy.^3.*uydx+(-1).*T.^3.*ux.*uy.^3.*uydx+2.*A.*P2.*uy.^4.*uydx+2.*P1.*uydy+2.*A.*P1.*uydy+(-1).*T.^3.*uydy+2.*P1.*ux.^2.*uydy+2.*A.*P1.*ux.^2.*uydy+(-1).*T.^3.*ux.^2.*uydy+(-2).*P2.*ux.*uy.*uydy+(-2).*A.*P2.*ux.*uy.*uydy+(-2).*P2.*ux.^3.*uy.*uydy+(-2).*A.*P2.*ux.^3.*uy.*uydy+2.*P1.*uy.^2.*uydy+2.*A.*P1.*uy.^2.*uydy+(-1).*T.^3.*uy.^2.*uydy+2.*P1.*ux.^2.*uy.^2.*uydy+2.*A.*P1.*ux.^2.*uy.^2.*uydy+T.^3.*ux.^2.*uy.^2.*uydy);

    R = [R1 R2 R3 R4 R5];


    [dt, fval, exitflag, output, jacobian] = fsolve(@eq, dt, optimset('Jacobian', 'on', 'DerivativeCheck', 'off', 'Display', 'off', 'Diagnostics', 'off'));

    DT  = dt(:,            1:  GRID_SIZE);
    Dux = dt(:,  GRID_SIZE+1:2*GRID_SIZE);
    Duy = dt(:,2*GRID_SIZE+1:3*GRID_SIZE);
    DP1 = dt(:,3*GRID_SIZE+1:4*GRID_SIZE);
    DP2 = dt(:,4*GRID_SIZE+1:5*GRID_SIZE);

    ut = sqrt(1 + ux.^2 + uy.^2);
    dtT  = (DT -ux.*Tdx -uy.*Tdy) ./ut;
    dtux = (Dux-ux.*uxdx-uy.*uxdy)./ut;
    dtuy = (Duy-ux.*uydx-uy.*uydy)./ut;
    dtP1 = (DP1-ux.*P1dx-uy.*P1dy)./ut;
    dtP2 = (DP2-ux.*P2dx-uy.*P2dy)./ut;


    T  = T  + Dt * dtT;
    ux = ux + Dt * dtux;
    uy = uy + Dt * dtuy;
    P1 = P1 + Dt * dtP1;
    P2 = P2 + Dt * dtP2;

    Tdx = dx(T);
    Tdy = dy(T);
    uxdx = dx(ux);
    uxdy = dy(ux);
    uydx = dx(uy);
    uydy = dy(uy);
    P1dx = dx(P1);
    P1dy = dy(P1);
    P2dx = dx(P2);
    P2dy = dy(P2);


    % dxdxux = dxdx(ux);
    % dxdyux = dx(dy(ux));
    % dydyux = dydy(ux);
    % dxdxuy = dxdx(uy);
    % dxdyuy = dx(dy(uy));
    % dydyuy = dydy(uy);
    % dxdxT = dxdx(T);
    % dxdyT = dx(dy(T));
    % dydyT = dydy(T);



    phydro();
    pause(0.003);
    disp(sprintf('time step: %04d,  %6.2fs , exit: %d, iterations: %d, funeval: %d', t, toc, exitflag, output.iterations, output.funcCount));
  end

















    phydro();
    %pause(0.033);

  function [res J] = eq(g)
    DT  = g(:,            1:  GRID_SIZE);
    Dux = g(:,  GRID_SIZE+1:2*GRID_SIZE);
    Duy = g(:,2*GRID_SIZE+1:3*GRID_SIZE);
    DP1 = g(:,3*GRID_SIZE+1:4*GRID_SIZE);
    DP2 = g(:,4*GRID_SIZE+1:5*GRID_SIZE);

    % dxdtT = dx(dtT);
    % dydtT = dy(dtT);
    % dxdtux = dx(dtux);
    % dydtux = dy(dtux);
    % dxdtuy = dx(dtuy);
    % dydtuy = dy(dtuy);

    L1 = (1/2).*(1+ux.^2+uy.^2).^(-1/2).*(6.*DT.*T.^2+8.*Dux.*P1.*ux+8.*Duy.*P2.*ux+6.*Dux.*T.^3.*ux+4.*DP1.*ux.^2+9.*DT.*T.^2.*ux.^2+(-8).*Duy.*P1.*uy+8.*Dux.*P2.*uy+6.*Duy.*T.^3.*uy+8.*DP2.*ux.*uy+(-4).*DP1.*uy.^2+9.*DT.*T.^2.*uy.^2);

    L2 = (1/2).*(1+ux.^2+uy.^2).^(-2).*(4.*Dux.*P1+4.*Duy.*P2+3.*Dux.*T.^3+4.*DP1.*ux+9.*DT.*T.^2.*ux+12.*Dux.*P1.*ux.^2+10.*Duy.*P2.*ux.^2+9.*Dux.*T.^3.*ux.^2+8.*DP1.*ux.^3+18.*DT.*T.^2.*ux.^3+8.*Dux.*P1.*ux.^4+6.*Duy.*P2.*ux.^4+6.*Dux.*T.^3.*ux.^4+4.*DP1.*ux.^5+9.*DT.*T.^2.*ux.^5+4.*DP2.*uy+(-4).*Duy.*P1.*ux.*uy+8.*Dux.*P2.*ux.*uy+3.*Duy.*T.^3.*ux.*uy+10.*DP2.*ux.^2.*uy+(-4).*Duy.*P1.*ux.^3.*uy+6.*Dux.*P2.*ux.^3.*uy+3.*Duy.*T.^3.*ux.^3.*uy+6.*DP2.*ux.^4.*uy+4.*Dux.*P1.*uy.^2+6.*Duy.*P2.*uy.^2+6.*Dux.*T.^3.*uy.^2+4.*DP1.*ux.*uy.^2+18.*DT.*T.^2.*ux.*uy.^2+12.*Dux.*P1.*ux.^2.*uy.^2+6.*Duy.*P2.*ux.^2.*uy.^2+9.*Dux.*T.^3.*ux.^2.*uy.^2+4.*DP1.*ux.^3.*uy.^2+18.*DT.*T.^2.*ux.^3.*uy.^2+6.*DP2.*uy.^3+10.*Dux.*P2.*ux.*uy.^3+3.*Duy.*T.^3.*ux.*uy.^3+8.*DP2.*ux.^2.*uy.^3+4.*Duy.*P2.*uy.^4+3.*Dux.*T.^3.*uy.^4+9.*DT.*T.^2.*ux.*uy.^4+2.*DP2.*uy.^5);

    L3 = (1/2).*(1+ux.^2+uy.^2).^(-2).*((-4).*Duy.*P1+4.*Dux.*P2+3.*Duy.*T.^3+4.*DP2.*ux+(-4).*Duy.*P1.*ux.^2+6.*Dux.*P2.*ux.^2+6.*Duy.*T.^3.*ux.^2+6.*DP2.*ux.^3+4.*Dux.*P2.*ux.^4+3.*Duy.*T.^3.*ux.^4+2.*DP2.*ux.^5+(-4).*DP1.*uy+9.*DT.*T.^2.*uy+4.*Dux.*P1.*ux.*uy+8.*Duy.*P2.*ux.*uy+3.*Dux.*T.^3.*ux.*uy+(-4).*DP1.*ux.^2.*uy+18.*DT.*T.^2.*ux.^2.*uy+10.*Duy.*P2.*ux.^3.*uy+3.*Dux.*T.^3.*ux.^3.*uy+9.*DT.*T.^2.*ux.^4.*uy+(-12).*Duy.*P1.*uy.^2+10.*Dux.*P2.*uy.^2+9.*Duy.*T.^3.*uy.^2+10.*DP2.*ux.*uy.^2+(-12).*Duy.*P1.*ux.^2.*uy.^2+6.*Dux.*P2.*ux.^2.*uy.^2+9.*Duy.*T.^3.*ux.^2.*uy.^2+8.*DP2.*ux.^3.*uy.^2+(-8).*DP1.*uy.^3+18.*DT.*T.^2.*uy.^3+4.*Dux.*P1.*ux.*uy.^3+6.*Duy.*P2.*ux.*uy.^3+3.*Dux.*T.^3.*ux.*uy.^3+(-4).*DP1.*ux.^2.*uy.^3+18.*DT.*T.^2.*ux.^2.*uy.^3+(-8).*Duy.*P1.*uy.^4+6.*Dux.*P2.*uy.^4+6.*Duy.*T.^3.*uy.^4+6.*DP2.*ux.*uy.^4+(-4).*DP1.*uy.^5+9.*DT.*T.^2.*uy.^5);

    L4 = (1/2).*T.^(-2).*(1+ux.^2+uy.^2).^(-1).*((-8).*DT.*P1+(-8).*A.*DT.*P1+4.*DP1.*T+4.*A.*DP1.*T+2.*Dux.*P1.*T.*ux+2.*A.*Dux.*P1.*T.*ux+(-4).*Duy.*P2.*T.*ux+Dux.*T.^4.*ux+(-16).*DT.*P1.*ux.^2+(-16).*A.*DT.*P1.*ux.^2+8.*DP1.*T.*ux.^2+8.*A.*DP1.*T.*ux.^2+2.*Dux.*P1.*T.*ux.^3+2.*A.*Dux.*P1.*T.*ux.^3+(-8).*Duy.*P2.*T.*ux.^3+(-2).*A.*Duy.*P2.*T.*ux.^3+Dux.*T.^4.*ux.^3+(-8).*DT.*P1.*ux.^4+(-8).*A.*DT.*P1.*ux.^4+4.*DP1.*T.*ux.^4+4.*A.*DP1.*T.*ux.^4+(-4).*Duy.*P2.*T.*ux.^5+(-2).*A.*Duy.*P2.*T.*ux.^5+2.*Duy.*P1.*T.*uy+2.*A.*Duy.*P1.*T.*uy+4.*Dux.*P2.*T.*uy+(-1).*Duy.*T.^4.*uy+(-8).*DT.*P2.*ux.*uy+(-8).*A.*DT.*P2.*ux.*uy+4.*DP2.*T.*ux.*uy+4.*A.*DP2.*T.*ux.*uy+10.*Duy.*P1.*T.*ux.^2.*uy+6.*A.*Duy.*P1.*T.*ux.^2.*uy+6.*Dux.*P2.*T.*ux.^2.*uy+(-1).*Duy.*T.^4.*ux.^2.*uy+(-8).*DT.*P2.*ux.^3.*uy+(-8).*A.*DT.*P2.*ux.^3.*uy+4.*DP2.*T.*ux.^3.*uy+4.*A.*DP2.*T.*ux.^3.*uy+8.*Duy.*P1.*T.*ux.^4.*uy+4.*A.*Duy.*P1.*T.*ux.^4.*uy+4.*Dux.*P2.*T.*ux.^4.*uy+2.*A.*Dux.*P2.*T.*ux.^4.*uy+(-8).*DT.*P1.*uy.^2+(-8).*A.*DT.*P1.*uy.^2+4.*DP1.*T.*uy.^2+4.*A.*DP1.*T.*uy.^2+4.*A.*Dux.*P1.*T.*ux.*uy.^2+2.*Duy.*P2.*T.*ux.*uy.^2+4.*A.*Duy.*P2.*T.*ux.*uy.^2+2.*Dux.*T.^4.*ux.*uy.^2+(-8).*DT.*P1.*ux.^2.*uy.^2+(-8).*A.*DT.*P1.*ux.^2.*uy.^2+4.*DP1.*T.*ux.^2.*uy.^2+4.*A.*DP1.*T.*ux.^2.*uy.^2+(-8).*Dux.*P1.*T.*ux.^3.*uy.^2+(-4).*A.*Dux.*P1.*T.*ux.^3.*uy.^2+4.*Duy.*P2.*T.*ux.^3.*uy.^2+2.*A.*Duy.*P2.*T.*ux.^3.*uy.^2+4.*Dux.*P2.*T.*uy.^3+2.*A.*Dux.*P2.*T.*uy.^3+(-8).*DT.*P2.*ux.*uy.^3+(-8).*A.*DT.*P2.*ux.*uy.^3+4.*DP2.*T.*ux.*uy.^3+4.*A.*DP2.*T.*ux.*uy.^3+(-4).*Dux.*P2.*T.*ux.^2.*uy.^3+(-2).*A.*Dux.*P2.*T.*ux.^2.*uy.^3);

    L5 = (1/2).*T.^(-2).*(1+ux.^2+uy.^2).^(-1).*(8.*DT.*P1+8.*A.*DT.*P1+(-4).*DP1.*T+(-4).*A.*DP1.*T+(-2).*Dux.*P1.*T.*ux+(-2).*A.*Dux.*P1.*T.*ux+4.*Duy.*P2.*T.*ux+(-1).*Dux.*T.^4.*ux+8.*DT.*P1.*ux.^2+8.*A.*DT.*P1.*ux.^2+(-4).*DP1.*T.*ux.^2+(-4).*A.*DP1.*T.*ux.^2+4.*Duy.*P2.*T.*ux.^3+2.*A.*Duy.*P2.*T.*ux.^3+(-2).*Duy.*P1.*T.*uy+(-2).*A.*Duy.*P1.*T.*uy+(-4).*Dux.*P2.*T.*uy+Duy.*T.^4.*uy+(-8).*DT.*P2.*ux.*uy+(-8).*A.*DT.*P2.*ux.*uy+4.*DP2.*T.*ux.*uy+4.*A.*DP2.*T.*ux.*uy+(-4).*A.*Duy.*P1.*T.*ux.^2.*uy+2.*Dux.*P2.*T.*ux.^2.*uy+4.*A.*Dux.*P2.*T.*ux.^2.*uy+2.*Duy.*T.^4.*ux.^2.*uy+(-8).*DT.*P2.*ux.^3.*uy+(-8).*A.*DT.*P2.*ux.^3.*uy+4.*DP2.*T.*ux.^3.*uy+4.*A.*DP2.*T.*ux.^3.*uy+16.*DT.*P1.*uy.^2+16.*A.*DT.*P1.*uy.^2+(-8).*DP1.*T.*uy.^2+(-8).*A.*DP1.*T.*uy.^2+(-10).*Dux.*P1.*T.*ux.*uy.^2+(-6).*A.*Dux.*P1.*T.*ux.*uy.^2+6.*Duy.*P2.*T.*ux.*uy.^2+(-1).*Dux.*T.^4.*ux.*uy.^2+8.*DT.*P1.*ux.^2.*uy.^2+8.*A.*DT.*P1.*ux.^2.*uy.^2+(-4).*DP1.*T.*ux.^2.*uy.^2+(-4).*A.*DP1.*T.*ux.^2.*uy.^2+(-4).*Duy.*P2.*T.*ux.^3.*uy.^2+(-2).*A.*Duy.*P2.*T.*ux.^3.*uy.^2+(-2).*Duy.*P1.*T.*uy.^3+(-2).*A.*Duy.*P1.*T.*uy.^3+(-8).*Dux.*P2.*T.*uy.^3+(-2).*A.*Dux.*P2.*T.*uy.^3+Duy.*T.^4.*uy.^3+(-8).*DT.*P2.*ux.*uy.^3+(-8).*A.*DT.*P2.*ux.*uy.^3+4.*DP2.*T.*ux.*uy.^3+4.*A.*DP2.*T.*ux.*uy.^3+8.*Duy.*P1.*T.*ux.^2.*uy.^3+4.*A.*Duy.*P1.*T.*ux.^2.*uy.^3+4.*Dux.*P2.*T.*ux.^2.*uy.^3+2.*A.*Dux.*P2.*T.*ux.^2.*uy.^3+8.*DT.*P1.*uy.^4+8.*A.*DT.*P1.*uy.^4+(-4).*DP1.*T.*uy.^4+(-4).*A.*DP1.*T.*uy.^4+(-8).*Dux.*P1.*T.*ux.*uy.^4+(-4).*A.*Dux.*P1.*T.*ux.*uy.^4+4.*Duy.*P2.*T.*ux.*uy.^4+2.*A.*Duy.*P2.*T.*ux.*uy.^4+(-4).*Dux.*P2.*T.*uy.^5+(-2).*A.*Dux.*P2.*T.*uy.^5);




    res = [L1 L2 L3 L4 L5] + R;

    % compute the gradients. hyhyhy

    dTL1 = (3/2).*T.^2.*(1+ux.^2+uy.^2).^(-1/2).*(2+3.*ux.^2+3.*uy.^2);
    dTL2 = (9/2).*T.^2.*ux;
    dTL3 = (9/2).*T.^2.*uy;
    dTL4 = (-4).*(1+A).*T.^(-2).*(P1+      P1.*ux.^2+P2.*ux.*uy);
    dTL5 =   4 .*(1+A).*T.^(-2).*(P1+(-1).*P2.*ux.*uy+P1.*uy.^2);

    duxL1 = (4.*P1.*ux+3.*T.^3.*ux+4.*P2.*uy).*(1+ux.^2+uy.^2).^(-1/2);
    duxL2 = (1/2).*(1+ux.^2+uy.^2).^(-2).*(4.*P1+3.*T.^3+12.*P1.*ux.^2+9.*T.^3.*ux.^2+8.*P1.*ux.^4+6.*T.^3.*ux.^4+8.*P2.*ux.*uy+6.*P2.*ux.^3.*uy+4.*P1.*uy.^2+6.*T.^3.*uy.^2+12.*P1.*ux.^2.*uy.^2+9.*T.^3.*ux.^2.*uy.^2+10.*P2.*ux.*uy.^3+3.*T.^3.*uy.^4);
    duxL3 = (1/2).*(1+ux.^2+uy.^2).^(-2).*(4.*P2+6.*P2.*ux.^2+4.*P2.*ux.^4+4.*P1.*ux.*uy+3.*T.^3.*ux.*uy+3.*T.^3.*ux.^3.*uy+10.*P2.*uy.^2+6.*P2.*ux.^2.*uy.^2+4.*P1.*ux.*uy.^3+3.*T.^3.*ux.*uy.^3+6.*P2.*uy.^4);
    duxL4 = (1/2).*T.^(-1).*(1+ux.^2+uy.^2).^(-1).*(2.*P1.*ux+2.*A.*P1.*ux+T.^3.*ux+2.*P1.*ux.^3+2.*A.*P1.*ux.^3+T.^3.*ux.^3+4.*P2.*uy+6.*P2.*ux.^2.*uy+4.*P2.*ux.^4.*uy+2.*A.*P2.*ux.^4.*uy+4.*A.*P1.*ux.*uy.^2+2.*T.^3.*ux.*uy.^2+(-8).*P1.*ux.^3.*uy.^2+(-4).*A.*P1.*ux.^3.*uy.^2+4.*P2.*uy.^3+2.*A.*P2.*uy.^3+(-4).*P2.*ux.^2.*uy.^3+(-2).*A.*P2.*ux.^2.*uy.^3);
    duxL5 = (-1/2).*T.^(-1).*(1+ux.^2+uy.^2).^(-1).*(2.*P1.*ux+2.*A.*P1.*ux+T.^3.*ux+4.*P2.*uy+(-2).*P2.*ux.^2.*uy+(-4).*A.*P2.*ux.^2.*uy+10.*P1.*ux.*uy.^2+6.*A.*P1.*ux.*uy.^2+T.^3.*ux.*uy.^2+8.*P2.*uy.^3+2.*A.*P2.*uy.^3+(-4).*P2.*ux.^2.*uy.^3+(-2).*A.*P2.*ux.^2.*uy.^3+8.*P1.*ux.*uy.^4+4.*A.*P1.*ux.*uy.^4+4.*P2.*uy.^5+2.*A.*P2.*uy.^5);

    duyL1 = (4.*P2.*ux+(-4).*P1.*uy+3.*T.^3.*uy).*(1+ux.^2+uy.^2).^(-1/2);
    duyL2 = (1/2).*(1+ux.^2+uy.^2).^(-2).*(4.*P2+10.*P2.*ux.^2+6.*P2.*ux.^4+(-4).*P1.*ux.*uy+3.*T.^3.*ux.*uy+(-4).*P1.*ux.^3.*uy+3.*T.^3.*ux.^3.*uy+6.*P2.*uy.^2+6.*P2.*ux.^2.*uy.^2+3.*T.^3.*ux.*uy.^3+4.*P2.*uy.^4);
    duyL3 = (1/2).*(1+ux.^2+uy.^2).^(-2).*((-4).*P1+3.*T.^3+(-4).*P1.*ux.^2+6.*T.^3.*ux.^2+3.*T.^3.*ux.^4+8.*P2.*ux.*uy+10.*P2.*ux.^3.*uy+(-12).*P1.*uy.^2+9.*T.^3.*uy.^2+(-12).*P1.*ux.^2.*uy.^2+9.*T.^3.*ux.^2.*uy.^2+6.*P2.*ux.*uy.^3+(-8).*P1.*uy.^4+6.*T.^3.*uy.^4);
    duyL4 = (-1/2).*T.^(-1).*(1+ux.^2+uy.^2).^(-1).*(4.*P2.*ux+8.*P2.*ux.^3+2.*A.*P2.*ux.^3+4.*P2.*ux.^5+2.*A.*P2.*ux.^5+(-2).*P1.*uy+(-2).*A.*P1.*uy+T.^3.*uy+(-10).*P1.*ux.^2.*uy+(-6).*A.*P1.*ux.^2.*uy+T.^3.*ux.^2.*uy+(-8).*P1.*ux.^4.*uy+(-4).*A.*P1.*ux.^4.*uy+(-2).*P2.*ux.*uy.^2+(-4).*A.*P2.*ux.*uy.^2+(-4).*P2.*ux.^3.*uy.^2+(-2).*A.*P2.*ux.^3.*uy.^2);
    duyL5 = (1/2).*T.^(-1).*(1+ux.^2+uy.^2).^(-1).*(4.*P2.*ux+4.*P2.*ux.^3+2.*A.*P2.*ux.^3+(-2).*P1.*uy+(-2).*A.*P1.*uy+T.^3.*uy+(-4).*A.*P1.*ux.^2.*uy+2.*T.^3.*ux.^2.*uy+6.*P2.*ux.*uy.^2+(-4).*P2.*ux.^3.*uy.^2+(-2).*A.*P2.*ux.^3.*uy.^2+(-2).*P1.*uy.^3+(-2).*A.*P1.*uy.^3+T.^3.*uy.^3+8.*P1.*ux.^2.*uy.^3+4.*A.*P1.*ux.^2.*uy.^3+4.*P2.*ux.*uy.^4+2.*A.*P2.*ux.*uy.^4);

    dP1L1 = 2.*(ux+(-1).*uy).*(ux+uy).*(1+ux.^2+uy.^2).^(-1/2);
    dP1L2 = 2.*ux.*(1+ux.^2).*(1+ux.^2+uy.^2).^(-1);
    dP1L3 = (-2).*uy.*(1+uy.^2).*(1+ux.^2+uy.^2).^(-1);
    dP1L4 = 2.*(1+A).*T.^(-1).*(1+ux.^2);
    dP1L5 = (-2).*(1+A).*T.^(-1).*(1+uy.^2);

    dP2L1 = 4.*ux.*uy.*(1+ux.^2+uy.^2).^(-1/2);
    dP2L2 = uy.*(1+ux.^2+uy.^2).^(-1).*(2+3.*ux.^2+uy.^2);
    dP2L3 = ux.*(1+ux.^2+uy.^2).^(-1).*(2+ux.^2+3.*uy.^2);
    dP2L4 = 2.*(1+A).*T.^(-1).*ux.*uy;
    dP2L5 = 2.*(1+A).*T.^(-1).*ux.*uy;



    J = zeros(3*GRID_SIZE*GRID_SIZE, 3*GRID_SIZE*GRID_SIZE);

    for i=1:GRID_SIZE*GRID_SIZE
      du = sparse(zeros(GRID_SIZE, GRID_SIZE));
      du(i) = 1;
      %dxdu = dx(du);
      %dydu = dy(du);

      Jxy = dTL1 .* du;
      J(                      1:  GRID_SIZE*GRID_SIZE,i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dTL2 .* du;
      J(  GRID_SIZE*GRID_SIZE+1:2*GRID_SIZE*GRID_SIZE,i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dTL3 .* du;
      J(2*GRID_SIZE*GRID_SIZE+1:3*GRID_SIZE*GRID_SIZE,i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dTL4 .* du;
      J(3*GRID_SIZE*GRID_SIZE+1:4*GRID_SIZE*GRID_SIZE,i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dTL5 .* du;
      J(4*GRID_SIZE*GRID_SIZE+1:5*GRID_SIZE*GRID_SIZE,i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);

      Jxy = duxL1 .* du;
      J(                      1:  GRID_SIZE*GRID_SIZE,1*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = duxL2 .* du;
	    J(  GRID_SIZE*GRID_SIZE+1:2*GRID_SIZE*GRID_SIZE,1*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = duxL3 .* du;
	    J(2*GRID_SIZE*GRID_SIZE+1:3*GRID_SIZE*GRID_SIZE,1*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = duxL4 .* du;
	    J(3*GRID_SIZE*GRID_SIZE+1:4*GRID_SIZE*GRID_SIZE,1*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = duxL5 .* du;
	    J(4*GRID_SIZE*GRID_SIZE+1:5*GRID_SIZE*GRID_SIZE,1*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);


      Jxy = duyL1 .* du;
      J(                      1:  GRID_SIZE*GRID_SIZE,2*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = duyL2 .* du;
	    J(  GRID_SIZE*GRID_SIZE+1:2*GRID_SIZE*GRID_SIZE,2*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = duyL3 .* du;
	    J(2*GRID_SIZE*GRID_SIZE+1:3*GRID_SIZE*GRID_SIZE,2*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = duyL4 .* du;
	    J(3*GRID_SIZE*GRID_SIZE+1:4*GRID_SIZE*GRID_SIZE,2*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = duyL5 .* du;
	    J(4*GRID_SIZE*GRID_SIZE+1:5*GRID_SIZE*GRID_SIZE,2*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);



      Jxy = dP1L1 .* du;
      J(                      1:  GRID_SIZE*GRID_SIZE,3*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dP1L2 .* du;
	    J(  GRID_SIZE*GRID_SIZE+1:2*GRID_SIZE*GRID_SIZE,3*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dP1L3 .* du;
	    J(2*GRID_SIZE*GRID_SIZE+1:3*GRID_SIZE*GRID_SIZE,3*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dP1L4 .* du;
	    J(3*GRID_SIZE*GRID_SIZE+1:4*GRID_SIZE*GRID_SIZE,3*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dP1L5 .* du;
	    J(4*GRID_SIZE*GRID_SIZE+1:5*GRID_SIZE*GRID_SIZE,3*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);

      Jxy = dP2L1 .* du;
      J(                      1:  GRID_SIZE*GRID_SIZE,4*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dP2L2 .* du;
	    J(  GRID_SIZE*GRID_SIZE+1:2*GRID_SIZE*GRID_SIZE,4*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dP2L3 .* du;
	    J(2*GRID_SIZE*GRID_SIZE+1:3*GRID_SIZE*GRID_SIZE,4*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dP2L4 .* du;
	    J(3*GRID_SIZE*GRID_SIZE+1:4*GRID_SIZE*GRID_SIZE,4*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dP2L5 .* du;
	    J(4*GRID_SIZE*GRID_SIZE+1:5*GRID_SIZE*GRID_SIZE,4*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
    end
    % sqrt(sum(sum(res.*res)))
  end


  function r = phydro()
    frame_active = (mod(t-1,int16(1/Dt/SECONDS_PER_LIGHTYEAR/30)) == 0) || EVERY_FRAME;
    if VIDEO & ~frame_active 
      return
    end
    if VIDEO
      disp(t);
    end
    plotsx = 4;
    plotsy = 3;

    subplot(plotsy, plotsx, 1);
    image((T-TEMP_OFFSET) * 128 + 128); 
    
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vt,frame);
    % end

    subplot(plotsy, plotsx, 2);
    image(ux*128 + 128);
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vvx,frame);
    % end

    subplot(plotsy, plotsx, 3);
    image(uy*128 + 128);
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vvy,frame);
    % end

    subplot(plotsy, plotsx, 4);
    image(sqrt(ux.*ux + uy.*uy)*128 + 128);
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vv,frame);
    % end

    subplot(plotsy, plotsx, 5);
    image(dtT+128);

    subplot(plotsy, plotsx, 6);
    image(dtux*4+128);
    subplot(plotsy, plotsx, 7);
    image(dtuy*4+128);

    subplot(plotsy, plotsx, 9);
    image(P1*16+128);
    subplot(plotsy, plotsx, 10);
    image(P2*16+128);

    %truesize(fig);
    if (VIDEO) & (frame_active)
      frame = getframe(gcf);
      writeVideo(vt,frame);
    end

    
    % subplot(plotsy, plotsx, 9);
    % image(dtdtvx*8 + 128);
    % subplot(plotsy, plotsx, 10);
    % image(dtdtvy*8 + 128);
    % subplot(plotsy, plotsx, 11);
    % image(Pi*2 + 128);
    % subplot(plotsy, plotsx, 12);
    % image(B*2 + 128);
    return;

  end

  % [-1 0 1
  function r = dx2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for y = 1:GRID_SIZE
      r(1,y) = - arg(GRID_SIZE,y)/2 + arg(2, y)/2;
    end
    for y = 1:GRID_SIZE
      r(GRID_SIZE,y) =  -arg(GRID_SIZE-1,y)/2 + arg(1, y)/2;
    end
    for x = 2:GRID_SIZE-1
      for y = 1:GRID_SIZE
        r(x, y) = (arg(x+1,y) - arg(x-1,y))/2;   
      end
    end 
    r = r / DX;
  end
  function r = dy2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for x = 1:GRID_SIZE
      r(x,1) = - arg(x,GRID_SIZE)/2 + arg(x,2)/2;
    end
    for x = 1:GRID_SIZE
      r(x,GRID_SIZE) =  -arg(x,GRID_SIZE-1)/2 + arg(x,1)/2;
    end
    for x = 1:GRID_SIZE
      for y = 2:GRID_SIZE-1
        r(x, y) = (arg(x,y+1) - arg(x,y-1))/2;   
      end
    end 
    r = r / DX;
  end
  function r = dx(arg)
    r = d1 * arg;
  end
  function r = dy(arg)
    r = arg * (d1.');
  end
  function r = dxdx(arg)
    r = d2 * arg;
  end
  function r = dydy(arg)
    r = arg * (d2.');
  end
  function r = dxdx2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for y = 1:GRID_SIZE
      r(1,y) = arg(GRID_SIZE,y) + arg(2, y) - 2*arg(1,y);
    end
    for y = 1:GRID_SIZE
      r(GRID_SIZE,y) =  arg(GRID_SIZE-1,y) + arg(1, y) - 2*arg(GRID_SIZE,y);
    end
    for x = 2:GRID_SIZE-1
      for y = 1:GRID_SIZE
        r(x, y) = arg(x+1,y) + arg(x-1,y) - 2*arg(x,y);
      end
    end 
    r = r / DX / DX;
  end
  function r = dydy2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for x = 1:GRID_SIZE
      r(x,1) = arg(x,GRID_SIZE) + arg(x,2) - 2*arg(x,1);
    end
    for x = 1:GRID_SIZE
      r(x,GRID_SIZE) =  arg(x,GRID_SIZE-1) + arg(x,1) - 2*arg(x,GRID_SIZE);
    end
    for x = 1:GRID_SIZE
      for y = 2:GRID_SIZE-1
        r(x, y) = arg(x,y+1) + arg(x,y-1) - 2*arg(x,y);   
      end
    end 
    r = r / DX / DX;
  end

end


